;-----------------------------------------------------------------
; 
; bivariate PDF distribution of RN222, SO2, and T  
;
;   - Using a user-defined function for cleaner code
;   - Subselecting an array using its coordinate arrays
;   - Using the "where" function
;   - Using "raster" contour mode for better delination of data
;   - Using "plt_pdfxy"
;-----------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
load "~/wrk/lib/lib_radon.ncl" 


begin 

  datapath = "./" 
  plotpath = "./" 
  plotform = "eps"     

  dist = "RN222_SO2" 
  dist = "RN222_T_P" 
  dist = "RN222_T" 

  sitename = "Freiburg" 
  sitename = "Beijing"
  sitename = "Schauinsland" 
  sitename = "Socorro" 
  sitename = "EUR" 
  sitename = "ASIA2" 
  sitename = "USA" 
  sitename = "AUS" 
  sitename = "PANNEL" 

  plotname = dist + "_" + sitename 

  flnma = "PDF_ASIA2.nc" 
  flnmb = "PDF_EUR.nc"
  flnmc = "PDF_USA.nc" 
  flnmo = "PDF_pannel.nc" 

  StringFontHeightF  = 0.035
  cnLineLabelFontHeightF = 0.025
  tmXBLabelFontHeightF = 0.032
  tmYLLabelFontHeightF = 0.032
  tiMainFontHeightF  = 0.05
  tiXAxisFontHeightF = 0.030 ;;5
  tiYAxisFontHeightF = 0.030 ;;5
  lbLabelFontHeightF = 0.030 ;;30
  gsnStringFont = "helvetica"

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/)

  abcd = abcd + ") "


  fla = addfile(datapath+flnma, "r")
  flb = addfile(datapath+flnmb, "r")
  flc = addfile(datapath+flnmc, "r")

  pdfaa = fla->PDFA
  pdfab = fla->PDFB
  pdfac = fla->PDFC

  pdfba = flb->PDFA
  pdfbb = flb->PDFB
  pdfbc = flb->PDFC

  pdfca = flc->PDFA
  pdfcb = flc->PDFB
  pdfcc = flc->PDFC


  system("rm "+flnmo) 

  flo = addfile(flnmo,"c") 

  flo->PDFAA=pdfaa 
  flo->PDFAB=pdfab 
  flo->PDFAC=pdfac

  flo->PDFBA=pdfba 
  flo->PDFBB=pdfbb 
  flo->PDFBC=pdfbc

  flo->PDFCA=pdfca 
  flo->PDFCB=pdfcb 
  flo->PDFCC=pdfcc

;-----------------------------------------------------------------
; plots
;-----------------------------------------------------------------

  colormap = "testcmap" 
  colormap = "amwg" 

  plot = new(9,"graphic")

  wks = gsn_open_wks(plotform, plotname)

  gsn_define_colormap(wks,colormap)

  res                      = True
  res@gsnDraw              = False    ; do not draw picture
  res@gsnFrame             = False    ; do not advance frame
;;----------------------------------------------------------------------
;; contour setting 
;;----------------------------------------------------------------------
  res@cnFillOn             = True     ; turn on color fill
  res@cnFillMode           = "RasterFill"       ; Raster Mode
  res@cnLinesOn            = False    ; no contour lines
  res@cnLineLabelsOn       = False    ; no contour line labels
  res@cnInfoLabelOn        = False
  res@gsnSpreadColors      = True     ; use full colormap
 ;res@cnLabelBarEndStyle   = "ExcludeOuterBoxes"

;;----------------------------------------------------------------------
;; labelbar 
;;----------------------------------------------------------------------
  res@lbLabelBarOn         = True       ; no individual label bar
  res@lbLabelAutoStride    = True
  res@lbLabelStride        = 1          ; every other label bar label
 ;res@lbOrientation        = "Vertical" ; vertical label bar
  res@lbLabelFontHeightF   = lbLabelFontHeightF
 ;res@pmLabelBarWidthF     = 0.5        ; default is shorter
 ;res@pmLabelBarHeightF    = 0.1        ; default is taller

;;----------------------------------------------------------------------
;; axis
;;----------------------------------------------------------------------
  res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) "
  res@tiXAxisString        = ""
  res@tiXAxisOffsetYF      = 0.01 
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF

;;----------------------------------------------------------------------
;; string
;;----------------------------------------------------------------------
  res@gsnLeftStringFontHeightF   = StringFontHeightF
  res@gsnCenterStringFontHeightF = StringFontHeightF
 ;res@gsnRightStringFontHeightF  = StringFontHeightF
  res@gsnStringFont = gsnStringFont

;;----------------------------------------------------------------------
;; tickmark
;;----------------------------------------------------------------------
  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; contour levels 
;;----------------------------------------------------------------------
  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  0.2               ; set min contour level
  res@cnMaxLevelValF       =  2.                ; set max contour level
  res@cnLevelSpacingF      =  0.2               ; set contour spacing
  ;res@cnMaxLevelValF       =  1.                ; set max contour level
  ;res@cnLevelSpacingF      =  0.1               ; set contour spacing

  ip = 0 

  res@gsnLeftString        = abcd(ip)+"ASIA    WCRP1995 " ;;""
  res@gsnCenterString      = "" ;;abcd(ip)+"WCRP1995 "
 ;res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = "Temperature (K)" 
  plot(ip) = gsn_csm_contour (wks, pdfaa, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"ASIA    SW1998 / 1.6 " ;;""
  res@gsnCenterString      = "" ;;abcd(ip)+"SW1998 / 1.6 "
 ;res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = ""
  plot(ip) = gsn_csm_contour (wks, pdfab, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"ASIA    MERGE " ;;""
  res@gsnCenterString      = "" ;;abcd(ip)+"MERGE " 
 ;res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = ""
  plot(ip) = gsn_csm_contour (wks, pdfac, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"EUR    WCRP1995 " ;;""
  res@gsnCenterString      = "" ;;abcd(ip)+"WCRP1995 "
 ;res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = "Temperature (K)"
  plot(ip) = gsn_csm_contour (wks, pdfba, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"EUR    SW1998 / 1.6 " ;;""
  res@gsnCenterString      = "" ;;abcd(ip)+"SW1998 / 1.6 "
 ;res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = ""
  plot(ip) = gsn_csm_contour (wks, pdfbb, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"EUR    MERGE " ;;""
  res@gsnCenterString      = "" ;;abcd(ip)+"MERGE "
 ;res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = ""
  plot(ip) = gsn_csm_contour (wks, pdfbc, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"NA    WCRP1995 " ;;""
  res@gsnCenterString      = "" ;; abcd(ip)+"WCRP1995 "
  res@tiYAxisString        = "Temperature (K)"
  res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  plot(ip) = gsn_csm_contour (wks, pdfca, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"NA    SW1998 / 1.6 " ;;""
  res@gsnCenterString      = "" ;; abcd(ip)+"SW1998 / 1.6 "
  res@tiYAxisString        = ""
  res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  plot(ip) = gsn_csm_contour (wks, pdfcb, res)

  ip = ip + 1 

  res@gsnLeftString        = abcd(ip)+"NA    MERGE " ;;""
  res@gsnCenterString      = "" ;; abcd(ip)+"MERGE "
  res@tiYAxisString        = ""
  res@tiXAxisString        = "IPRR (cm~S~-2~N~ s~S~-1~N~ ) " 
  plot(ip) = gsn_csm_contour (wks, pdfcc, res)


  resP                     = True                ; modify the panel plot
  ;resP@gsnMaximize         = True        ; must include w/ Paper Orientation
  resP@gsnPanelYWhiteSpacePercent = 0 ;;5
  resP@gsnPanelXWhiteSpacePercent = 0 ;;5
  resP@gsnPanelBottom      = 0.05 
 ;resP@gsnPanelTop         = 1.00 
 ;resP@gsnPanelLeft        = 1.00 
 ;resP@gsnPanelRight       = 0.00 
 ;resP@txString            = title
 ;resP@gsnPanelRowSpec     = True                   ; tell panel what order to plt
 ;resP@gsnPanelLabelBar    = True
 ;resP@lbLabelAutoStride   = True
 ;resP@lbOrientation       = "vertical"          ; vertical label bar


  gsn_panel(wks,plot,(/3,3/),resP)

;-----------------------------------------------------------------
; Use "plt_pdfxy" to plot
;-----------------------------------------------------------------

;  res@gsnCenterString      = "PDF"
;  plot(0) = plt_pdfxy (wks, pdfa, res)
;
;  resP                     = True 
;  resP@txString            = "" 
;  resP@gsnPanelRowSpec     = True              
;  gsn_panel(wks,plot,(/1,1/),resP)

  system("ps2pdf "+plotname+".eps") 

end

                 
